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Abstract 

An efficient time-stepping algorithm is proposed based on operator-splitting and the space-time discontinuous 
Galerkin finite element method for problems in the non-classical theory of thermoelasticity. The non- 
classical theory incorporates three models; the classical theory based on Fourier’s law of heat conduction 
resulting in a hyperbolic-parabolic coupled system, a non-classical theory of a fully hyperbolic extension, 
and a combination of the two. The general problem is split into two contractive sub-problems, namely the 
mechanical phase and the thermal phase. Each sub-problem is discretised using space-time discontinuous 
Galerkin finite element method resulting each to be stable which then leads to unconditional stability of the 
global product algorithm. A number of numerical examples are presented to demonstrate the performance 
and capability of the method. 

Keywords: Operator-splitting, Space-time discontinuous Galerkin finite element, Non-classical theory of 
thermoelasticity, Fourier’s law, Second sound, Contractive. 


1. Introduction 

In some solids thermal energy can be transmitted by the mechanism of wave-like propagation of heat, 
unlike the usual mechanism of conduction by diffusion. This phenomenon of heat conduction as waves, 
known as second sound , has been observed experimentally (see, for example, mm for an extensive survey 
of experimental works involving propagation of heat as a thermal wave). 

The classical theory of heat conduction based on Fourier’s law, fails to model the second sound phe¬ 
nomenon. Moreover, the classical theory permits infinite speed of propagation of parts of a localized initial 
heat pulse, which is paradoxical from a physical point of view. As a result, efforts have been made in 
an attempt to find a consistent model of heat propagation that is capable of capturing the second sound 
phenomenon with finite speed (see, for example, mm for a review of models of heat conduction as waves). 

One of the alternative theory for formulating the propagation of heat in a general way that aims at 
capturing the second sound phenomenon was proposed by Green and Naglidi mmmm. The theory of 
Green and Naghdi is based on three types of constitutive relations for the heat flux: Type I is equivalent 
to the classical theory based on Fourier’s law. Type II permits the propagation of a localized heat signal as 
thermal wave without dissipation (see [5] for a remark on the appropriateness of this classification). Type 
III is the most general theory, which includes both type I and II as special cases, in which second sound 
phenomenon is supported while dissipation is incorporated in the process. 

The thermomechanical coupling of non-classical heat conduction with classical elasticity is the subject of 
non-classical thermoelasticity. Extensive overviews of the non-classical thermoelasticity of Green and Naghdi 
can be found in pa inim- Theoretical results concerning the non-classical theory have been addressed in 
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several research works. In |12] exact solutions are obtained for thermal wave propagation in one dimension. 
Results on the existence and uniqueness of solutions of non-classical problem of thermoelasticity can be 
found, for example, in [13| and the references therein. 

Designing a robust and efficient numerical solution strategy for strongly coupled problems of hyperbolic- 
type is challenging. This is particularly the case for the non-classical theory of thermoelasticity where 
hyperbolic (or nearly hyperbolic) heat conduction equation is coupled with the classical hyperbolic elasticity 
problem. A standard approach for solving such time-dependent problems is the Method of Lines (MoL) in 
which the governing partial differential equation is first discretised in space using the finite element method 
(FEM) leading to a system of ordinary differential equations, which can then be solved using the finite 
difference method. Despite its popularity, MoL struggles to accurately solve problems involving propagation 
of sharp gradients or discontinuities mm- 

Recently, a great deal of attention has been invested in designing a spatial Discontinuous Galerkin (DG) 
approach for convection-dominated problems; see for example [16]. However, these methods, like MoL, 
are based on decoupling space and time in the sense that space and time are treated differently. Hulbert 
and Hughes mug introduced a powerful scheme based on a space-time DG finite element methodology 
for linear elastodynamics problems. In their approach, space and time are treated simultaneously and the 
unknown fields are allowed to be discontinuous in time while continuous in space. Recently, the space-time 
DG method has been used in m for classical thermoelasticity, using a monolithic approach where all the 
unknown fields are solved for simultaneously. 

Recently, in 1181 a numerical solution approach based on MoL was proposed for non-classical thermoelas¬ 
ticity in which time integration was done in two ways: continuous Galerkin FEM for type II and III, while 
mixed-discountinuous Galerkin FEM for the classical problem based on the Fourier’s law of heat conduction. 
In their approach a streamline-upwind numerical stabilization was added to localize numerical oscillations 
due to the propagation of sharp thermal wave. 

In the current work, we extend the existence and uniqueness results otained in m for type II theory to 
the more general problem of type III thermoelasticity. We also present a novel numerical algorithm for the 
non-classical thermoelasticity based on an operator-splitting technique motivated by Armero and Simo m 
for classical thermoelasticity, coupled with a space-time DG methodology that extends the work of Hulbert 
and Hughes [14] which was formulated for linear elastodynamics. The major contributions of this work are 
twofold: (i) the adaptation of the operator-splitting strategy for classical thermoelasticity first proposed by 
m to the non-standard theory, in which the operator defining non-standard thermoelasticity is split in a 
way that the resulting sub-operators retain the same contractive behaviour as the global operator; and (ii) 
the time-DG formulation in which continuity of the unknown fields is enforced weakly by using an L 2 -inner 
product in contrast to the energy-norm used in |14j . 

The rest of this paper is organized as follows. In Section [2] the governing equations of the non-classical 
theory are summarized in a general framework of type III thermoelasticity. Well-posedness and physically 
meaningful boundary and initial conditions are also discussed in this section. An operator-splitting strategy 
for the problem of type III thermoelasticity is proposed and the resulting sub-operators are analysed in 
Section]!] In Section]!] time-DG formulation is proposed for the sub-problems and stability of the individual 
algorithms and the global one is analysed in detail. A number of numerical examples both in 1 D and 2-D 
are presented in Section[5]to demonstrate the excellent performance and capability of the proposed numerical 
scheme. Finally, concluding remarks and some open problems are discussed in Section [6] 

2. Model problem: Non-classical thermoelasticity (type III) 

This section summarizes the equations governing the non-classical theory of thermoelasticty of type III 
as proposed by Green and Naghdi. Well-posedness of the problem is analysed. Results obtained here will 
serve as the basis of design and analysis of the numerical algorithm that will be presented in later sections. 

Governing equations 

Let H C K d , with 1 < d < 3 be the reference placement of a continuum body B with smooth boundary 
r. Following Green and Naghdi’s theory of thermoelasticity of type III, the system of partial differential 


2 


equations governing the thermomechanical interaction in the solid B are 


u = v \ 

pi) = div[Ce(ti) — mi?] + pb 

infix!, (1) 

a = 0 

pci) = div[k 2 Va + k 3 V0] — 0om : e{u) + pr. 

where I = [0, T] is the time interval of interest of length T > 0. Superimposed dots denote time derivatives. 
The displacement and the velocity vector fields are denoted by u and v respectively. The scalar field i? 
denotes the relative temperature with respect to a uniform reference value 0 O > 0 such that the absolute 
temperature 0 is given by 0 = i) + 0o- The quantities b and r are the prescribed body force and heat 
source. 

Green and Naghdi’s theory of non-classical thermoelasticity is based on the inclusion of a state variable, 
known as the thermal displacement a, that is defined in terms of an empirical temperature T (which is 
assumed to coincide with the absolute temperature 0) through equation ®3, see, for example, [20j H3 [21] 
and the references therein. 

The symbol e(u) = sym(Vit) denotes the small strain tensor associated with a displacement u. It is 
assumed that the elasticity tensor C has the following properties: 

C-ijkl = Cjikl = Cijik, ( 2 ) 

Cijkl = Cfcirp (3) 

^ijkitijCki > 0 for any non-zero symmetric tensor e. (4) 

Equations |2| and ^ are minor and major symmetries of C, while equation Q is the positive definiteness 
of C. The coupling second-order tensor m is of the form 


m = 3 wk1, 


where u, k = A + 2/3/r, and 1 denote, respectively, the thermal expansion coefficient, the bulk modulus and 
the identity second-order tensor, and p and A are the Lame constants. It is assumed that the tensor k 2 is 
symmetric and positive-definite, and that k 3 are symmetric and positive-semidefmite. The scalars p > 0 
and c > 0 denote the material density and heat capacity. 

Remarks: 


1. The non-classical theory of thermoelasticity of type III ([I]) is the most general one in that it contains 
both type I and II as a special cases. If k 2 Va is omitted form ([!]), then one obtains type I (or classical 
thermoelastic model) where a parabolic heat conduction equation is coupled with the hyperbolic 
mechanical equation. On the other hand, if k 3 is set to zero, one obtains the type II thermoelastcity 
where, now, part of the system ([T]), that is responsible for heat conduction, is hyperbolic (non-classical 
heat conduction). 

2. Under the assumption of mechanical and thermal isotropy the elasticity tensor C, and the tensors k 2 
and k 3 become 

C = A1 <g) 1 + 2 //I, k 2 = fc 2 l, and k 3 = fc 3 1, 
where I denotes the forth-order identity tensor, and /c 2 > 0, and fc 3 > 0 are constants. 


3. The free energy i(>, and hence the stress tensor er and the entropy density ?? are given by 

pi/j = : Ce - dm : £ - ^ - i)S 0 , 

Z Z t/Q 

T = ^ = c (u) - mtf, and m = + m :*(«) + S„, 

ae o0 0o 


( 5 ) 


where Sq is the absolute entropy density. 
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4. The heat flux vector q within the non-classical theory of thermoelasticity of type III is defined as 

q = — [k 2 Va + k 3 V0]. 

Using the entropy constitutive relation ©3 the coupled system ([I]) can be written in terms of q as 

u = v 

pi) = div[Ce — mil] + ph 
a = 0 

pQoV = div[k 2 Va + k 3 V0] + pr j 

It is this form of the dynamical system which is crucial in designing the computational scheme based on 
operator-splitting in latter sections. 


infix [0, T 1. 
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2.1. Initial and boundary conditions 

Let {r u ,T t } and {r^,T g } be two partitions of T, each contains mutually disjoint subsets; that is, 

r = r u ur t = r tf ur„ with r u n r t = n r, = 0. 

Let u : T u xI-> R d , t : T t x I —> R d , d : T# x I —> R, and q : T q x I —> R be prescribed displacement, 

traction, thermal displacement and flux fields. Thus the boundary conditions are given by 

u = u on T u x I, d = d on T# x I, 

(7) 

ern = t on T t x I, q ■ n = q on T ? x I, 

where n denotes the outward unit normal field to T. It is easy to observe the analogy between the two set of 

equations: the mechanical part Qi ,2 and the thermal part |l]) 3j 4 . In such analogy, we clearly see that the 
displacement u goes with the thermal displacement a (in fact, it is this analogy that motivated the name 
thermal displacement 0), and the velocity v goes with the absolute temperature 0, and hence with d. As 
a consequence, however, one would expect a thermal Dirichlet boundary condition is given in terms of a 
as it is customarily the case in the mechanical part that u is prescribed as Dirchlet boundary condition. 
The thermal Dirichlet boundary condition, in this case, is given via the relative temperature d (and hence 
absolute temperature 0). The reason for this is that, usually, boundary conditions are prescribed in terms 
of physical quantities, which can be measured, which, in the thermal case, is the relative temperature, d (or 
0 ). ' 

Furthermore, the initial conditions read 

u(x, 0) = m°(x), v(x, 0) = v°(x), 
a(x, 0) = a°(x), d(x, 0) = d°(x), 

where u°, v°, a 0 , and d° are prescribed initial displacement, velocity, thermal displacement and absolute 
temperature respectively. In prescribing an initial thermal state, a thermal configuration is assumed so that 
the initial thermal displacement a is homogeneous, that is a 0 = 0, while, the physically observable quantity, 
the relative temperature can be initiated at a non-zero value. 


2.2. Well-posedness: Dissipation and conservation 

Let L c , T c , M Cl and K c be characteristic scalar quantities with the dimensions of length, time, mass, 
and temperature, respectively. Define the dimensionless variables as 



' 1' 
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After introducing the dimensionless variables, the non-dimensional form of (jT]) become 

u = v, 

pv = div[Ce(ii) — mi9] + ph , 
a = 0, 

pcQ = div[k 2 Vb + k 3 V0] — ©orh : e(u) + pf , 


where the spatial and time derivatives are with respect to the dimensionless space and time variables, and 

n2l 


C = 
k 2 = 


\L c Tl 1 

C, 

m = 

\L c TlKA 


\T 2 1 
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If we drop the bars in the notations of equation Q, similar expressions as in equation ([I]) is obtained. For 
the reminder of this section, whenever the system 0 is mentioned, unless stated otherwise, it refers to its 
non-dimensional form, and the initial and boundary conditions should also be understood accordingly. 

The positive-definiteness property of C and k 2 , and the positive-semidefiniteness of k 3 imply that the 
system (jT|) together with the initial and boundary conditions ([7]) and <[8j) define an evolution equation of a 
general form 

X(t) = Ax(t) +f] 

X(0) = X° I 


in V, 


( 10 ) 


where A is a closed linear operator with dense domain D( A) C V defined in some suitable Banach space 
V. For the case of non-classical linear thermoelasticity, for the sake of simplicity, we consider homogeneous 
Dirichlet boundary condition with respect to both u and a and the space V 


V:= {(m, v, a, 0) T e [H 1 (fi)] d x [L 2 (fi)] d x H^fi) x L 2 (M) : u = 0, a = 0 on T} , 


(ID 


is a Hilbert space. 

The abstract solution vector x = ( M , v , a, 0) T £ V, while the linear operator A and the source term f 
in (10) are defined by 


A*:= 


-div[Ce(u) — null 
P 

0 

—div[k 2 Va + k 3 V01 — —m : s(v) 
.pc pc 



' 0 ' 


b 

, f:= 

0 


1 


-r 

Lc -I 


( 12 ) 


We consider an inner product, (•, -) v on V defined by 

(X>x)v = ( e ( u )> Ce(w)) + (pv, v) + (k 2 Va, Vd) + d) . 


(13) 


where (■, •) denotes the standard L 2 -inner product pairing of tensor, vector, or scalar fields, that should be 
understood in context, and k^ = k 2 pc/0o and c* = pc/0o- The norm on V induced by the inner product 
") v is denoted by || • || V - 

Note that the linear differential operator A : V(A) C V —> V is closed and the space 


[Hj(fi) n H 2 (fi)] d X [H),(H)] d x (^(O) n H*(rt)) x H'{n) c v(A), 


(14) 


is dense in V. Hence T>(A) is dense in V. 

A very important inequality concerning the evolution equation (101 is dissipativity property of the defining 
operator A. An operator A on closed subspace D( A) of a Hilbert space V endowed with an inner product 
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(-,-) v is said to be dissipative if it satisfies the inequality (Ax, x)v — f° r eac ^ X € T>( A) |]jjl- If the 
operator A is dissipative, the norm of the solution of the corresponding evolution equation is monotonically 
decreasing in time, which is referred to as contractivity of the solution. That is, for a solution x of the 


evolution problem (10), assuming dissipativity of A and f = 0, we have 


4llxl|v = 4 (x,x)v = 2 (x,x)v = 2 ( A X, x) < 0. 


dt 


dt 


(15) 


Now, we shall show that the operator A that defines the problem (10) is dissipative. Let x = ( u , v , a, id) 


be in the domain of A, D( A), satisfying the homogeneous boundary condition. Then 

(X, A x)v = (e(u),Ce(v)) + (pv, -div[Ce(it) - mi?) 

P 

+ (k*Va, V0) + (c*i?, —div[k 2 Va + k 3 V0] - —m : e(v)) 
pc c 

= (e(tt), Ce(v)} - (e(v), Ce(u)) + (e(»),ini?) + (^ L k 2 Va, V0) 

0o 

- (^k 2 V0, Vo) - (^k 3 V0, V0) - < 1 ?, m : e(v)> 

Wo Wo 

= -(^k 3 V0,V0) <0. 

Wo 


(16) 


In the general context, since k 3 is positive-semidefinite equation (16) leads to dissipation. In the limiting 


case where k 3 = 0 (type II) the above argument implies conservation of energy-norm 

S(t) := My = \ f [e(u) : C e(u) + pv ■ v + k(jVa • Va + c*^ 2 ]dfl. 
2 Jn 


(17) 


This is the reason why type II is also referred to as the theory of thermoelasticity without energy dissipation , 
see, for example, mmm- 

Another important relation concerning the operator A is that it should satisfy the following: for all 
X* € V, there exists x i n 2?(A) such that 

X-Ax = x\ ( 18 ) 

in other words, the operator (1 — A) : T>( A) —► V is onto. 


To show that A satisfies the relation (18), we proceed as follows: let x = (u,v,a,'d) T and %* = 
(it*, v*, a*, d*) T then from the definition of A equation (18), implies that 


u — V = u 


v -div[Ce(it) — mi?] = v*, 

a — i? = a *, 

i?-^-div[k 2 Va + k 3 V0] + —m : e(v) = i?*. 

pc c 


(19) 


Substitution of equations (|19|i and (19) 3 into the remaining equations of (19) leads to the (equilibrium) 


problem: find x = ( u , v , a, t?) J £ D(A) such that v = u — it*, i? = a — a* and satisfying 

p 2 0oit — p©odiv[Ce(it) — ma] = u, 
pea — div[kVo:] + pQ om : e(u) = a, 


( 20 ) 


where u = p 2 0 o it*+p 2 0oii*+p0odiv[mQ!*], a = pca*+pcd* — div[k 3 Va*]+p©om : e(u*), and k = k 2 +k 3 . 

The weak form of equation (20) reads: find x = (it, v, a, i?) T £ V such that v = u — u*,fl = a — a* and 
satisfying 


B( X ,o = m 

6 


(21) 










for all £ = (w, is, P, w) £ V. The bilinear form B (•, •) and the right hand side functional l(-) are given by 


B(x,£) = (p 2 Qou,w) + (pQ 0 Ce(u),e(w)) - (p0 o ma,£(ir)) 
+ (pea, p) + (kVa, V/3) + (p0 o m : e(u), p), 

KQ = (u,w) + (a, P). 


( 22 ) 

(23) 


in equation (23) represents duality 


Note that u £ [H -1 (f2)] d and a £ H _1 (n) and here the symbol ( 
pairing in their respective spaces. 

From the definition of B(-, •), we can easily see that it is a bounded bilinear form. Since 

B(x,X) = (p 2 Q 0 u,u) + (p<d 0 Ce(u),e(u)} + (pea,a) + (kVa, Va), 


(24) 


then B(-, •) is ([Hj(f2)]“ x HQ(fl))-elliptic. By applying Lax-Milgram theorem we conclude that there exists 
X £ V which solves the weak problem (21), and hence solves equation (19). Therefore, this proves the 


ontoness of the resolvent operator (1 — A). 

In conclusion, we have seen that the operator A defining the non-classical linear thermoelasticity (type 
III) 

i) is closed, 

ii) has dense domain T>( A) in V, 

iii) is dissipative, and 

iv) is such that (1 — A) : V(A) C V —» V is onto. 

Therefore, by the Lumer-Phillips theorem, A generates a strongly continuous semigroup of contractions 


in V, see, for example [13] and the references therein. In other words, the problem (10) is well-defined 


and contractive. This also means that the dynamical system represented by the equation of non-classical 
thermoelasticity of type III is, in general, stable in the sense of Lyapunov. 


3. Algorithms based on operator-splitting strategy 


as 


Consider an abstract evolutionary problem of the form ( 10 ). Assume that A can be expressed additively 

(25) 


A — Ai + A 2 , 

such that the operators Aj, i = 1, 2 define sub-problems 

Xi =AiX! X,(0)=X°, * = 1,2. 


(26) 


Let A ^ be consistent and stable time-stepping algorithms corresponding to the sub-problems (26) with At 
representing the time step length that the algorithms step up the state vectors, Xu from a given time t 
to t + At. A time-stepping algorithm, A At , for the global problem is obtained by taking products of the 
algorithms as 

(27) 


A At = A: 


At 


, At 


o Af' 

The algorithm A At is referred to as Lie- Trotter-Kato product formula |22j . It is sometimes called a sequential 
split algorithm or single pass algorithm. In addition to the discretisation error in the individual algorithms 
A At , the application of the operator-splitting strategy introduces another source of error known as the 


splitting error. The splitting error associated to the Lie-Trotter-Katto product formula (27) is of order of 
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magnitude O(At) (see, for example, [23]). It means that A At is only first order accurate. Higher-order 
algorithms based on operator-splitting strategies include: 


Marchuk-Strang split: 
Double-pass split: 


A a ‘ = Af t/2 o A At o Af' 5/2 
A a ‘ = i(A At o A At + A At o A At ). 


(28) 


The Marchuk-Strang split is second-order accurate, while the double-pass split is only first-order. Note that 
Lie-Trotter-Kato product formula depends on the order of operations, for example in (27), A At is applied 
first followed by A At , while the order of operations does not matter in the other two operator-splitting 
algorithms given in (28). 


3.1. Operator-splitting for non-classical thermoelasticity 

The non-classical theory of thermoelasticity is a coupling of two dynamical systems: the classical linear 
elasticity and the non-Fourier thermal conduction. A naive splitting of the system ([Tj) into a mechanical 
problem under constant thermal states (isothermal) and a thermal problem with a fixed configuration will 
result in at most a conditionally stable algorithm even if the sub-algorithms for the two processes are 
unconditionally stable [19:. Rather, care must be taken in splitting the two systems; this is usually dictated 
by an understanding of the underlying physics. In this respect, it makes sense if we split the operator A in 


(121 so that in the mechanical phase the entropy is held fixed (isentropic) while the temperature is allowed 


to vary, and in the thermal phase heat is allowed to be conducted while the configuration is fixed. In fact, 
in this split, it can be shown that each sub-process defines an evolution problem which is contractive, as is 
the global problem. Furthermore, consistent and stable algorithms for the sub-processes render a consistent 
and stable algorithm for the global problem by the way of operator-splitting strategy. 

To this end, inspired by the work of Arrnero and Simo [Tj3], taking the 4-tuple £ = (u 1 v,a,p) T as the 
state variables we consider the splitting of the system of equation ([!]) into 


u = v, 

pi) = div[Ce(w) — mi?] + pb, 
a = 0, 
pQar) = 0 , 


and 


u = 0 , 
pi) = 0, 
a = 0 , 

pQoV = div[k 2 Va + k 3 V0] + pr. 


(29) 


This corresponds to the additive splitting of the operator A = Ai + A 2 in (10) 


A,S = 


-div[Ce(u) — mi?] 

0 

0 


A 2 £ = 


0 

0 

0 

— divfeVa 
L Pc 


k 3 V0] 


(30) 


Using the same calculation as for the dissipativity of A in (16), we obtain the estimates such that for each 

£, 

<A 1 S,S) v = 0, 


(A 2 S, E)v = -(-^k 3 V0, V0) < 0. 

w 0 


(31) 


In general case, since k 3 is positive-semidefinite, both operators A 3 and A 2 are dissipative. In particular, if 


k 3 = 0, the systems that accounts for thermal conduction, (29)2, is energy conserving, i.e. it represents heat 


conduction without energy loss in a rigid body. Moreover, it can be shown that both of the sub-operators 
satisfy additional conditions in order to generate strongly continuous semigroups of contraction just like how 
it was done in Section [2721 

Now, having two sub-operators Ai and A 2 generating contractive semigroups let us assume that there 
corresponds two algorithms A At , *=1,2 which are B-stable ( non-linearly stable ) ; that is 


l|AfV-Af t xll<llx n -All 


for all x n , X e V. 


(32) 












In the linear case this means that each of the algorithms satisfy the estimate (see [23 and the references 
therein): 

||A^||v<l, i = 1,2. (33) 

where V* is the dual of V. Let {x n }„ 6 N and {x"} ne N two sequences in V generated by the Lie-Trotter- 
Kato product formula A At corresponding to two initial conditions x(0) = x° and x(0) = X° respectively. 
Then the product formula A At satisfies the stability estimate 

iix n+1 - r +1 iiv = iia a v - A At riiv 

= \\A At [A At x n }-A At [A At x n \\v 

<||AfV-Af‘x"|| v (||Af|| v <l) 

= llx n -X n ||v (||Af*||v < 1), (34) 

which proves the non-linear stability of the global algorithm corresponding to the product formula A At . 

The numerical scheme that is going to be formulated in the subsequent sections is based on the operator¬ 
splitting approach and time-discontinuous Galerkin finite element method. A Lie-Trotter-Kato product 
formula is applied to merge algorithms for the two phases. Hence, the product formula can also be viewed 
in the sense of predictor-corrector regime, where the sub-algorithm for the mechanical phase is used as 
predictor and that of the thermal phase as a corrector. 

4. Time-discontinuous Galerkin finite element method (T-DG FEM) 

Let Th = {fi e } be a triangulation of H, where H denotes the closure (the union of the interior and 
boundary) of O, such that 

f2= \J n e . (35) 

finely 

Denote the space of scalar piecewise polynomials on the mesh Th by £? 3 h \ 

= Wn e c°m ■■ <p h \ a . G P j (n% !! e er4, (36) 

where P- 7 '(fi e ) denotes the set of polynomials of degree at most j defined on f l e . 

Consider a partition of the time domain I = [0, T] into the collection {!„ = [t n , t n +i]} n ~Q, N £ N of 
non-overlapping subintervals. The time step length is = t n+ \ — t n for n = 0,1, • • • , N — 1 with 

0 = to < ti < • • • < = T. 

For each time sub-domain I n we consider the space-time domain of the form Q n = Cl x I n referred to as the 
n th space-time slab. 

Admissible scalar functions, (j ) h , that we consider in T-DG FEM will be polynomials in time t £ I n with 
coefficients from the spatial function space ^ >3 h \ i.e. 

<t> h {x,t) = Y^${x)t l , (37) 

i 

where t l is a monomial in t £ I n of order i £ N. Denote by S h (Q n \ j , /), the space of admissible functions on 
the space-time domain Q n of degree j + l (that is, j in space and l in time): 

S h {Q n \j,l ) = jV : <p h {x,t ) = ^2<p^(x)t\ $ £ &° h , (x,t) £ Q„| . (38) 

In fact, it is easy to observe that the space S h (Q n ;j,l) is generated by the tensor products of the basis 
elements of the spaces £? 3 h and P l {I n )~ the set of polynomials in time of degree at most j. 
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Remark: 

1. The space-time mesh for Q n is composed of cells with one element thickness in the time direction; 
and in each cell in the slab, time and space are orthogonal to each other, i.e. each cell is of the 
form Tt e x I n . Nevertheless, the formulation being developed here can be easily modified in terms of 
non-orthogonal space-time elements and with slabs composed of more than one element thickness in 
time direction as well. 

2. The approach would readily accommodate the use of an adaptive mesh refinement procedures. In 
such cases, there may be cells with time direction thickness less than A t n embedded in slab Q n . In 
this case, at each hanging node the solution must be constrained so that the hanging nodes will be 
condensed out later. 


Notations 

Let tp and ip be functions defined on space time domain slab Q n . Some frequently used notations are 

a) Spatial L 2 inner product 

(V?,V>):= / dO. 

J n 

b) Space-time L 2 inner product 

J (‘£,'0) d t. 

c) Space-time boundary integrals on = T t x /„ and F n = T q x I n 


= [ f dT d t, 

■> r rl J r t 

(ip,ip) Fri = J J (pip dT d t. 


c) Right/left limit of a discontinuous function in time at t n 

(fi{t±):= lim ip(t n +e). 

e->0± 

d) Temporal jump of a discontinuous function at t n 

Hn = - <p{tn)- 

4-1. Mechanical problem 

In the mechanical phase the entropy is held fixed, so that the last equation of the left hand system in 


(291; that is the equation 


P 0°?7 = ^ [pcd + pOom : e(u) + 5 0 ] = 0 on VL x [t n , t n+ 1 ) 


(39) 


is solved in closed form to obtain an intermediate temperature d 1 . This leads to an explicit formula for d 1 
in terms of the state variables at time step t n from the left and at time value t £ (t n , t n + 1 ), that is, 


d J (f) = d(t n ) - —m : e{u(t) - u(t n )) for t G (t n ,t n+1 ). 
c 


(40) 


We substitute this result into the mechanical problem (29) i to obtain 

u = v, 

pi) = div[C ad £(u)] + /, 
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( 41 ) 




where C ad = C + (0o/c)m ® m, in the terminology used in m, is referred to as adiabatic elasticity tensor 
and / = pb — m[x}(t~) + (0o/c)m : u(t~)] and d(t~) and u{t~) denote the temperature at the end of the 
previous space-time slab, Q n -i- Note that the adiabatic elasticity tensor C ad remains positive-definite and 
symmetric as C. 

On the current space-time slab, Q n , we use same boundary conditions as given in 0 for the mechanical 
fields but the initial conditions, in this case, are the solution for u and v at the end of the previous slab; 
that is u{t~ ) and v(t~). 

To define the T-DG FEM formulation of the mechanical problem (411, we first define the trial and weight 
function spaces for displacement u and velocity v vector fields as 


T d = {u h £ [S h {Q n -,j, l)] d -u h =u on T u X /„} , 
T" = { v h £ [S h (Q n ;j,l)] d : v h = u onT u x /„} , 
w” = {w h £ [S h (Q n ;j,l)] d : w h = 0 on T u x /„} , 
K = {<p h € [S h (Q n \j, l)] d :<p h =0 on T u x I n ) , 


where l~ h and T h , W f and W are trial and weight function spaces for displacement and velocity vector 
fields respectively. The T-DG FEM is formulated as: find U h = ( u h ,v h ) T £ T h x T h such that for all 
V h = (■ w h ,ip h ) T £ W“ x Wl 

A"(U\V h ) = b"(V h ), (43) 


where 


4 , (u\v h ) 


{u h , w h ) Qn - ( v h , w h ) Qn + (, pv h , <p h ) Qn + (C ad E(u h ), e{ip)) Q „ 
+ ( uh (tl), w h {tl)) + (pv h {t+), v h {tl)) 


KiV 11 ) = + (/> l fi h )Qn + {u h {t n ),w h {t+)) + (pv fl (t n ) 1 ip' l (tl)). 


■, h (++\ 


(44) 


Remark: 


(1) The main difference between the DG formulation presented here in (43) and that of Q3] is the inner 
product used to enforce the equation of motion (41) weakly. In our formulation we use the L 2 -inner 


product to weakly enforce the mechanical problem while in na an energy-inner product is used. 


(2) The formulation (43) is consistent in the sense of a time-stepping algorithm. This can be seen from 
the Euler-Lagrange form of (|43| given by 


0 = AZ(U h ,V h )-b“{V h ) 

= ( u h — v h , w h )Q n + ( pv h — div[C a(i £(^^) ,^ ] — /, ip h ) (equation of motion) 

+ ([it ,l ] rM w h (tD) (displacement continuity) 

+ (lpv h ]n, V h {ti)), (velocity continuity) (45) 


that upon substitution of a sufficiently smooth solution pair ( u,v) T of the strong form (41) into (45), 
the weak forms of the jumps and the equation of motion vanish. 


(3) The jump terms are used to improve the stability of the scheme without degrading the accuracy. As 
a result, the formalism can be readily extended to the non-linear case without eliminating the jump 
term from the displacement-velocity relation. 

(4) One of the consequences of using the L 2 -inner product is that a Dirchlet-type boundary condition may 
not be necessary to define the velocity trial and weight function spaces. Instead we can use 


r h =wi = [s h (Q n -,j,i)r 


(46) 
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4-1.1. Stability: The mechanical algorithm 

For the sake of simplicity, we assume a homogeneous source term / = 0 and boundary conditions, i.e. 
t = 0 and u = 0. We claim that the formulation (431 renders an unconditionally stable time-stepping 
algorithm. That is to say: 

Vn = 0,1, • • • ,N —1, 

where $ M (U(t)) is the total mechanical energy of U = (u,v) T at time t, given by 


(47) 

(48) 


(49) 


^m( C 7 W) = [e(«W) ■ C ad e{u(t)) + pv(t) ■ v(t)\ da 

For the analysis we use elliptic and L 2 interpolation operators tv : [ifr„(^)] d — > kV h (t) and ft : [L 2 (fl)] d 
W“(t) respectively defined as: for u £ [F7 ru (f2)] d and w £ [L 2 (Q)] d 

{C ad e{Tvu)],£(<p h )) = (C ad e{u),E{ip h )), V<P h e W“(i), 

(ftru,r/^) = (u>,i/> h ), £ W“(t), 

where W"(i) referees to the space of functions in W” at a fixed but arbitrary t £ We also use the fact 
that [ 25 

div[C Qd e(7rM)] = frdiv[C ad e(«)]. (50) 

Now, given the solution U h (t4) = (u h (t~),v h (t~)) T of (43) at the end of the previous space-time slab, 
Qn-i, and let U h (t ) = (u h (t), v h (t)) T be the solution of (43) in the current space-time slab, Q n . Replace 
V h = (ftdiv[C ad e(M ?! ')], 0) T in (43) to obtain 

0 = (u h ,ivdiv[C ad £{u h )])Q n - (v h ,Tvdiv[C ad £{u h )}) Qrb 

+ (I« ?l ]n,7tcliv[C ad £(M ,l (t+))]), 


(51) 


Applying (50) in (51) and the definition of the projection operator tv and using integration by parts (note 
the homogeneous boundary conditions) leads to 

0 = (e(u h ), C ad e{u h )) Qn - {e{v h ), C ad e(u h )) Qri 


+ ([£(«'*)]», C od £(« h (i+))). 


Again substituting V h = ( 0,v h ) T into (43) yields 

0 = ( pv h , v h ) Qn + (C ad £(u h ),e{v)) Qn + {\pv h ]„, v h (t+)). 


(53) 


Adding the equations ( |52| and (53) we obtain 

(e(u h ),C ad e(u h )) Qn + ( pv h ,v h ) Qri + (le(u h )] n ,C a . d e{u h {t+))) + ([ pv h j n ,v h {t+)) = 0. (54) 


Taking the time derivative out of the space integral, (541 becomes 

^(£(“' l (^+i)).C od e(M' t (f” +1 )))Q„ + ^(pv h {t~ +1 ),v h (t~ + i))q„ 

- \^{u h (t+)),C ad £{u h (t+))) Qn - ^(pv h {t+),v h (t+)) Qri 

+ (le(u h )UC ad e(u h (t+))) + (lpv h j n ,v h (tt)) = 0. 

After some algebraic manipulation we obtain 

4(^(t n p)) + 4(n) = s M {u h {t~)). 


(55) 


(56) 


Since <t> M {\U h \ n ) is non-negative the energy equation (56) leads to the estimate (47) that renders the 
scheme for the mechanical phase (43) unconditionally stable. In fact, the total numerical dissipation added 
is precisely equal to 


IV-l 




n =0 
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Remark: 


The use of projection operators in (49) and (51) reveals an important point, that is, the current T-DG 
formulation can be converted into the one given in [14j . 


4-2. Thermal problem 

The solution U h (t~ +1 ) of the mechanical phase is known at the end of the current space-time slab. The 
objective, in the present phase, is to solve for the thermal states II(t) = (a(t), , t G I n . Hence, the 

global solution at the end of the current slab will be (U h (t~ +1 ),Tl h (t~ +1 )). 

The operator-splitting is performed based on the state vector X = (it, v, a , p) T . As a result, we enforce 
the problem in the thermal phase using the conservation form 


a = 0 , 

p&oV = div[lt 2 Va + k^VG] + pr. 


(57) 


cp 

Here, recall that the entropy density p is obtained from the relation pp = + m : e(u) + Sq. 

G 0 

In the thermal phase, the displacement u and the velocity v are fixed at the corresponding values at the 
end of the current slab in the mechanical phase; that is, 


u T (t) = u M (t n+1 ), v T {t) =v M (t n+ 1 ), t€l n , 


(58) 


where the subscripts T and M represents the values of the fields in the thermal and mechanical phases, 
respectively. Thus, the time derivative of the terms m : e(u) + Sq vanishes, and consequently the left hand 
side of equation (57)2 becomes 


p@oP = pcd. 

In addition, the jump in the entropy, in this case, reads as 


(59) 


IpQovjn = P®ov(t+) - p<dov(t~) 

= + P® 0 m : s{u M (t~ + 1 )) + So] - [pcd{t~) + pG 0 m : e(u(t~)) + S 0 ] (60) 

= \pcti\n + P©om : [s(u M (t~ + 1)) - e(w(t~))]. 

It should not cause any confusion if we drop the superscript M in the equation ( |60| so that the jump term 
can be written as 

Ip® ovin = {pc&ln + P©om : [e(u(t~ +x )) - e(u(t~))]. (61) 

To define the T-DG formulation for the thermal phase we first define the thermal displacement and 

q, i9 a t9 

temperature trial and weight function spaces T h , T h and W h , W h respectively based on S h (Q n ' j,l) and 
the boundary condition requirements. 


7]= {a h e S h (Q n ;j,l) : a h = G on T Q x /„} , 
7f = G S h (Q n ;j,l) : G k = 6 on T a x I n } , 
= {p h G S h {Q n ]j, l) : 0 h = 0 on T a x I n } , 
W,;’ = {a h G S h (Q n ;j,l) :a h = 0onT a x /„} . 


Formally, the T-DG FEM formulation of the thermal phase on the domain Q. n is defined as: find n /l = 
(a, , d) T gTj” x such that for each A h = (/ 3 h , a h ) T G x 

Al(U h ,A h ) = b T ri (A h ), (63) 


13 






where 


< (n\ A h ) = (a\p h ) Qn - (e h , /3%„ + {pcd h , a h ) Qn + ([k 2 Va h + k 3 V9 h ],Va h ) Qn , 

+ (a\t+),P h (t+)) + (pc1) h (t+),a h (t+)) 

K(A h ) = {a h (t~),P h (t+)) + {pcd h (t~),o h (t+)) + (p© 0 [e(u(t“)) - e(u(t~ + 1 ))],a h (t~)) 
+ (h,a h ) Frl + ( pr , a h ) Qn . 


The relation between the DG formulation (631 and the point-wise form (57) is apparent from the Euler- 
Lagrange form 


0 = Al(n h ,A h ) - bjA h ) 

+ {a-<d h ,/3 h ) Qn 

+ (pQoi) + di v[k 2 Va h + k 3 Vf? ,i ] + pr, a h )Q n (Equation of Motion) (65) 

+ (la h ln,/3 h (t+)) (a-continuity) 

+ (bc^Jn, a h (t +)), (^-continuity) 


which reveals that a sufficiently smooth solution of the strong problem (57) also satisfies (63), and vice versa, 
while the jump terms are vanished at the smooth solution. This also proves the consistency of the T-DG 
scheme of the thermal problem. 

The unconditionally stability of the scheme (63) can also be shown along the same line of argument used 
for the mechanical case. 


Remark: 

Again, the use of the L 2 -inner product to enforce the thermal problem allows one to omit the boundary 
restriction when we define the thermal displacement trial and weight function space, i.e. 

X :=S h (Q n -,j,l)=:W a h . (66) 

This is a very important observation in terms of practical implementation. 

As we have seen from Section[3]that consistent and stable sub-algorithms render a consistent and stable global 
algorithm in the sense of time-stepping algorithms based on operator-splitting. Both the mechanical and the 
thermal phase algorithms are shown to be consistent and unconditionally stable. Therefore, the algorithm 
for the global problem based on Lie-Trotter-Kato product formula is consistent and unconditionally stable. 
Moreover, the convergence of the global scheme follows from the well known result stated below. 

Theorem 1 (Lax Equivalence Theorem). For consistent numerical approximations, stability is equiv¬ 
alent to convergence. 


5. Numerical results 

In this section, we present a range of results for type II and III problems of non-classical thermoelasticity. 
We start by comparing convergence of the proposed splitting scheme against a monolithic approach in which 
all the governing equations are discretised simultaneously using the time-DG finite element method. For 
this, a 1-D problem of non-dimensional form is considered. The result shows excellent agreement between 
the monolithic and the splitting scheme. Then we go on to present various results in 1-D and 2-D. The 
examples in this case are designed to illustrate two key features of the time-DG scheme: (i) its performance 
in solving problems that involves the propagation of sharp gradients without creating spurious oscillations; 
and (ii) its capability in capturing the unique aspects of non-classical theory, for example, propagation of 
thermal wave and its complex response due to the coupling of elasticity problem. 
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The family of problems considered in this section are organized as follows. To analyse the rate of con¬ 
vergence and capability of the proposed scheme, non-dimensional form of a 1-D non-classical thermoelastic 
problem is presented in Section |5.1[ The performance of the splitting algorithm is examined in Section [5.2| 
for an initial temperature pulse propagation in a two dimensional square plate under plane strain condition. 
Finally, in Section [5.3[ a quasi-static expansion of a thick walled, infinitely long cylinder in plane strain con¬ 
dition is presented, which is modelled as type I and type III theory of thermoelasticity and the remarkable 
difference of thermal responses between the two models are also analysed. 


5.1. Non-dimensional 1-D GNT 

The non-dimensional form of 1-D GNT problem given in ([!]) is 


d T u = v, 

d T v = -d\+b, 

d T a = d, 

d T '& = <9 ? [d^a + kd^fl] — e 2 d s v + s, 


with the dimensionless parameters 


e = (1L) 2 
1 ’ ’ 


0 O m 2 E 
pc 


and k 


h 

yfpe.' 


(67) 


( 68 ) 


where s 1 denotes the square of ratio of uncoupled velocities of the mechanical wave (or first sound) and 
thermal wave (or second sound), e 2 denotes the strength of the thermomechanical coupling, k represents 
the non-dimensional classical heat conductivity. The speed of first sound Vf is actually the speed of sound 
in the medium, that is 


v f = 



(69) 


where E denotes the Young’s modulus of the medium, while that of the second sound v s is a characteristic 
feature of the theory of non-classical heat conduction by Green and Naghdi that represents the speed in 
which a thermal disturbance travels through the medium: 


v s = 



(70) 


The non-dimensionless variables are given by 


£ = x 




T = t 1 t , U = U 1 M, 


V = - V, 

U 


a = a 1 a, 


■&= — ■&, 
a 


(71) 


where x c , t c , u a , a c are characteristic quantities having the same dimension as x, t, u, a respectively that 
can be chosen according to the relations 


x c u o mv f 

t c a c p 


(72) 


From the above equations we can observe that there are infinitely many ways of choosing the characteristics 
constants without changing the form of the system ( [67] ). 

The nondimensional energy counterpart of also referred to as the H 1 -norm, is given by 


f(x) = 



e 1 [d s u] 2 -\-v 2 + — [9 { d] 2 + y-'d 2 

Sr. Sr. 


d£, 


(73) 


and the L 2 -norm 



where x = (u, v, a , id) T is the state vector at a given time. 
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5.1.1. Convergence 

For the purpose of the convergence analysis an exact solution to problem ( |67| ) is obtained in such a way 
that source terms b and s are suitably prescribed such that a given state vector x = (u, d, d, d) T be an exact 
solution HZI- To this end let the source terms be 


b = — [(Ej — 1) sin(7r£) sin(7rr) + cos(7r£) cos(7tt)] , 


(75) 


s — —j— [/c7rsin(7r^) cos(7rf) + e 2 cos(7t£) cos(7rr)], 


so that the exact solutions are 


«(£, t) = a = - sin(7r£) sin(7rr), 

_ - 7T 

v(£, t ) = $ = sin(7r£) cos(ttt), 


(76) 


defined on the space-time domain (£, t) € [0, L\ x [0,T]. For convergence analysis the values of the non- 
dimensional parameters at taken at e 1 = 4, e 2 = 0.2, k = 0. Such set of values represents a strongly coupled 
problem of two purely hyperbolic systems (type II thermoelasticity). The space-time domain corresponds 
to L = 1 and T = 0.25. Bilinear finite element functions, Q 1, are used in each space-time slab with each 
element having an aspect ratio of one (i.e. h = = At). 

Fig. |T] (a) reports the spatial convergence results of the monolithic and operator-splitting approaches 
with error norms of the approximate solutions at t = T are computed using the H 1 - and L 2 -norms as 
given in equations (73) and (741, respectively . By construction the monolithic algorithm is only first-order 


accurate. Remarkably, it is shown that the error norms for both approaches are seen to overlap showing the 
error associated with the splitting is almost negligible. This demonstrably shows an increment in efficiency 
of the splitting algorithm, while maintaining the accuracy of the monolithic scheme. 

While Fig. [T] (b) presents the result of temporal convergence results of the two approaches with errors of 
the approximate solutions computed at the mid-point, £ = L/ 2, and r = T using the £2 vector norm. Again 
errors corresponding to both of the approaches are almost identical. Similarly, this shows the temporal error 
associated with the operator-split strategy is minimal. 


5.1.2. Laser pulse propagation 

Consider a one-dimensional bar occupying the interval £ € [0,1], heated by a pulsing laser applied at the 
left end having the form similar to the one considered in |26j for non-Fourier heat conduction problem: 



where D is the depth of the pulse, and r p is characteristic duration of the pulse. The bar is clamped at 
both ends at all times and with homogeneous initial conditions. We consider a situation in which a highly 
localized thermal pulse both in space and time described by the constants r p = 0.01 and D = 0.02 is 
applied at the left end of the bar. The parameters considered here are eq = 9, which represents 3 : 1 ratio 
of uncoupled speeds of first sound to second sound, and £2 = 1 accounting for a strongly coupled system. 
Bilinear elements are used in each space-time slab with mesh dimension A£ = At = h. The simulations are 
carried out over the period of T = 1 unit of non-dimensional time. The mesh parameter h = 0.001 is chosen 
such that the width of the pulse is greater than the mesh size. In other words, the mesh is chosen so that 
the thin laser pulse can be described accurately by the the bilinear finite elements. 

Fig. [ 2 ] (a) and (b) show the propagation, in space and time, of the thermal disturbance caused by the 
pulsing laser heat source applied at the left end of the bar, computed using the monolithic and the splitting 
schemes, respectively. As can be seen from the figures, immediately after the pulse is applied, two thermal 
waves with different amplitude and speed start to emerge. The bigger and the slower wave is the one which is 
driven by the thermal equations, while the smaller and the faster one is induced by the mechanical equations 
through the coupling. The bigger thermal wave travels with a speed slightly less than that of second sound; 
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(a) Spatial convergence 



(b) Temporal convergence 

Figure 1: Type II thermo-mechanical problem: Rate of convergence using monolithic and splitting ap¬ 
proaches where the error norms are computed at r = 0.25 over the whole spatial domain. 
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whereas, the smaller thermal wave is travelling with a speed slightly greater than that of first sound. For 
this reason, it appears that the larger wave traverses the bar once, while the smaller traverses it more than 
three times. Note that the ratio of uncoupled speed of first to second sound is exactly 3 : 1. 

Here there are two features which show the strength of the thermomechanical coupling: the first one 
is that the ratio of the speeds of the two thermal waves is noticeably different from what is expected in 
uncoupled case, and the other is that the coupling is strong enough to induce considerably large stress wave 
which in turn induce the faster thermal wave. 

Moreover, this problem represents a strongly coupled problem of two second-order hyperbolic problems 
involving propagation of sharp gradients. Such a problem is typically very difficult to approximate using 
the standard semi-discrete approach (MoL) unless some kind of stabilization term (or an artificial viscosity) 
is added , which is basically equivalent to changing the system from non-dissipative to dissipative, or very 
fine mesh is used together with very small time-step, which is undesirable from a computational cost point 
of view. 

What is remarkable about the current scheme is that it resolves the propagation of high gradients 
accurately while the amplitude of the thermal waves appear to be constant showing a very small numerical 
dissipation is added enough to damp out any numerical oscillation that could happen. The two approximate 
solution profiles Fig. [2] (a) and (b) are nearly identical. The agreement demonstrates that the splitting 
scheme maintains the accuracy of the monolithic scheme while the efficiency is considerably improved by 
the splitting scheme since two smaller systems are solved at each space-time slab. The result obtained here 
can be qualitatively compared to the one obtained in 1201 . 

As shown from Fig. [3j other than some small numerical instabilities when the waves interact either with 
the boundary or each other, the energy gained computed using the iJ 1 -norm, remains essentially constant 
after the pulse is applied. This phenomenon is characteristic feature of type II thermoelasticity which is 
proved in Section 2.2 While the L 2 -norm shows more profound variation than the energy-norm right after 
the pulse is applied and when the two waves interact each other but it shows no change when the waves 
interact with the boundary. These observation suggests that the numerical instability that is arisen from 
the the interaction of waves with the boundaries may come from errors in the gradient of the approximate 
solution states. 

Fig. [4] (a) and (b) show the temperature profiles for the same problem above but with k = 0.1, which cor¬ 
respond to Type III thermoelasticity, approximated using the monolithic and splitting schemes, respectively. 
This case is characterized by dissipation of energy while a wave scenario is still evident. The thermal wave 
driven by the temperature equations is damped out quickly, where as, the mechanically induced thermal 
wave remains localized for almost the entire duration and is travelling with speed nearly equal to the speed 
of the first sound. 


5.2. Two dimensional problem: Initial heat pulse propagation 


In this problem, we consider a non-dimensional form of type III problem of initial thermal pulse propa¬ 
gation in a square plate occupying the region Q = [—1,1] x [—1,1] under plane strain assumption. A similar 
problem with the dimensions is solved in [20]. The boundary of the specimen is mechanically clamped and 
fixed at the reference temperature @0 = 1 (i.e. the temperature of the ambient space). Initially, it was at 
rest but a temperature pulse is initialized at the center of the plate, i.e. the initial condition for the relative 
temperature 1 ) be 


d(x, 0) = A exp 


x ■ x 
D 


(78) 


where D , as in the previous example in Section |5.1.2[ is a constant characterizing the width of the initial 
temperature pulse and A is the amplitude. The material parameter used in the simulation are scaled 
according to the specifications summarized in Table [l] The time-DG finite element mesh consists of 8 node 
isoparametric cubes with one element thickness At = 0.01 in time direction and 100 x 100 spatial elements 
per each slab are used to sufficiently describe the initial thermal pulse propagation. 

Fig- E shows snapshots of propagation of an initial temperature pulse with D = 100 and A = 4 at times 
t = 0, t = 0.2, t = 0.3, and t = 0.4. The initial pulse may be thought of as a thermal configuration just after 
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£-axis 

(a) Monolithic 


m. 




£-axis 
(b) Splitting 

Figure 2: Propagation of laser pulse in type II thermoelasticity: temperature profile of the rod over the time 
period with e 1 = 9, e 2 = 0.5, k = 0 and = At = 0.001. 
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(a) Monolithic 



(b) Splitting 


Figure 3: Propagation of laser pulse in type II thermoelasticity: the ^-Energies and L 2 -norms correspond¬ 
ing to using monolithic and splitting approaches. 
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(b) Splitting 


Figure 4: Propagation of laser pulse in type III thermoelasticity: temperature profile of the rod over the 
time period with k = 0.1, e, = 9, e 2 = 0.5, and = At = 0.001 
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t = 0 


t = 0.2 



0.129 el 
0.127 el 
0.126 el 
0.124 el 
0.122 el 
0.121 el 
0.119 el 
0.117 el 
0.11 6 el 
0.114 el 
0.112 el 
0.111 el 


0 


t = 0.3 


t = 0.4 




0.118 el 
0.117 el 
0.11 6 el 
0.114 el 
0.113 el 
0.112 el 
0.111 el 
0.110 el 
0.108 el 
0.107 el 
0.106 el 
0.105 el 


0 


Figure 5: Temperature distribution in a square plate according to type III thermoelasticity where an initial 
pulse localized in space is initiated at the center. 
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Table 1: Initial pulse propagation: material properties 


Speed of first sound 

V E /p 

1.96 

Speed of second sound 

\J k 2 /p c 

0.65 

Conductivity ratio 

K 2 /K 3 

100 



Figure 6: Finite element mesh for the problem of expansion of a thick-walled cylinder. 


an intense and highly localized laser heat source is applied at the center. The temperature profile gradually 
widens and a smaller but faster mechanically driven wave emerges, while the the second sound wave is driven 
by the temperature equations moves with a slower speed. In this case, the classical conductivity parameter 
K 2 gives additional stability but it is not so high to smear out the two wave phenomena. 


5.3. Quasi-static case: Expansion of a thick walled cylinder 

This problem deals with the quasi-static thermo-mechanical interaction in a thick walled cylinder as it 
expands as a result of an inner wall Dirichlet-type boundary condition, in the linear and plane strain case. 
The material considered is isotropic both thermally and mechanically. The thermal variation is purely the 


result of mechanical changes (the expansion of the cylinder) unlike in the previous examples (Sections 5.1.2 
and 5.2) in which thermal variations cause mechanical effects. 

The cylinder has cross section occupying the region Q = {(x,y) : 7q < x 2 + y 2 < R 2 } with inner and 
outer radii ro = 10 mm and R = 20 mm, respectively. A zero heat flux boundary condition is maintained on 
the inner wall, while the outer wall is kept at the reference temperature 0o- The inner wall is dynamically 
prescribed a radial displacement of 1 mm per second, while the outer wall is mechanically free. 

The problem is analysed for 20 seconds until the inner wall reaches a radius of r = 3ro- The time- 
DG finite element mesh consists of trilinear shape functions of 56 elements around the circumference of 
the cylinder by 8 elements radially with one element thickness in the temporal direction with step length 
At = 0.1 s for each space-time slab. In this quasi-static case, since only the thermal equations contains 
temporal derivatives then the thermal fields are allowed to be discontinuous while the displacement field is 
continuous across the interfaces of each space-time slab. This implies that the numerical dissipation comes 
from the weak enforcement of the continuity of the thermal fields only. 

We consider two cases: the first is classical or type I thermoelasticity with = 45 N/sK and the other 

is type III thermoelasticity with k 2 = 90 N/K and = 30 N/sK. 

Fig- 0 shows temperature variations over time for each case sampled at the equally spaced points along 
the radial direction labeled A-E as shown in the Fig. [6] As expected, in both cases, the temperature of the 
entire cylinder is converging to the reference temperature as time increases. The sinusoidal thermal response 
of type III is due to a temperature wave moving back and forth indicating second sound phenomenon. 
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Figure 7: Temperature profile of five points in the cylinder which are 0 mm, 2.5 mm, 5 mm, 7.5 mm, and 
10 mm away from the inner wall and shown with the labels A, B, C, D, and E. 


6. Conclusion 

An operator-splitting strategy coupled with a space-time discontinuous Galerkin finite element method 
for the solution of transient and fully coupled initial-boundary problem of generalized thermoelasticity was 
presented. Well-posdeness of the problem in the general setting (type III) is proven using the theory of 
semigroups. The defining operator is split additively so that the first sub-operator represents an isentropic 
(adiabatic) elasticity in which the entropy density is held fixed, and the other is a non-standard heat 
conduction at fixed configuration. Both of the sub-problems are also shown to inherit the same contractivity 
property as the full problem. 

Each sub-problem is then discretised separately using a time-discontinuous Galerkin finite element 
method where the unknown fields are allowed to be discontinuous along the interfaces of each space-time 
slab. Weak continuity of the unknown fields is enforced using an L 2 -inner product which differs from the 
original time-discontinuous formulation using an energy-inner product fl5l rm which was formulated for 
linear elastodynamics problem. The unconditionally stability behaviour of each of the algorithms is proven 
without the need to add extra ‘artificial viscosity’. The algorithm for the global problem is finally obtained 
by way of Lie-Trotter-Kato product formula, leading to an unconditional stability. 

The results presented in this paper are demonstrated by a number of numerical examples in both one 
and two dimensional cases. The efficiency of the current numerical scheme were examined by comparing 
the rate of convergence of with the corresponding monolithic approach. The result shows that the splitting 
scheme not only it retains the accuracy of the monolithic scheme but also it improves the efficiency as 
two smaller problems are solved sequentially at each time-step. The capability of the splitting algorithm 
is tested using problems involving propagation of heat waves driven by a pulsing laser heat source and an 
initial temperature disturbance in one and two dimensions respectively. Furthermore, the capability of the 
non-standard thermoelasticity and the proposed numerical method to model the phenomenon of second 
sound in some solids is demonstrated by considering the quasi-static expansion of an infinitely long thick 
walled cylinder in plane stress. 
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The DG formulation proposed in this work may be extended to the non-linear regime without the need 
to eliminate the displacement-velocity relation in the formulation. Hence, a full recovery of the numerical 
dissipation in the non-linear case is possible. This will be the subject of a forthcoming work in which issues 
such as non-linear stability and the existence of Lyapunov function are discussed. 
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